Hannahglover/Glowworm
Mapping genes to single cell RNAseq

Package index
Search the Hannahglover/Glowworm package
Functions
9
Source code
14
Man pages
5
Home
/
GitHub
/
Hannahglover/Glowworm
/

In Hannahglover/Glowworm: Mapping genes to single cell RNAseq

Load packages, save & load data

library(Seurat)
library(ggpubr)
library(dplyr)
library(plyr)
library(tidyr)
library(harmony)
library(matrixStats)
library(patchwork)
library(DoubletFinder)
options(future.globals.maxSize = 8000 * 1024^2)

save(list=c("Islet.combined", "Islet.combined_UMAP", "GeneLists", "CheckUMAP", "DEGs", "ClusterFunc_All_RNA", "GenerateMetaData", "Islets_Barcs", "Beta_Clust",  "Acinar_Clust",  "Ductal_Clust",  "Alpha_Clust",  "Macro_Clust",  "Delta_Clust", "Gamma_Clust",  "Epsilon_Clust",  "Cycling_Clust", "WBC_Clust",  "Endo_Clust", "Stellate_Clust", "Islets_Barcs_List", "Beta_Barcs",  "Acinar_Barcs", "Ductal_Barcs",  "Delta_Barcs", "WBCs_Barcs","Stellate_Barcs", "Endo_Barcs", "Cycling_Barcs", "Epsilon_Barcs", "Gamma_Barcs", "Macro_Barcs", "Alpha_Barcs", "List_Seu", "Alpha_Seu",  "Beta_Seu",  "Delta_Seu", "Epsilon_Seu", "Gamma_Seu", "Acinar_Seu",  "Ductal_Seu", "Stellate_Seu",  "Endo_Seu",  "WBCs_Seu", "Macro_Seu", "List_Meta", "Cycling_Seu"), file = "~/Islet_StdIntegration_DbltRm.RData")
load("~/Library/CloudStorage/Box-Box/HG2553 Main Folder/Collins Lab Collab/Islet_StdIntegration_DbltRm.RData")

library(devtools)
install_github("campbio/celda", build_vignettes = TRUE)

vignette("decontX")

library(celda)
library(SingleCellExperiment)

sce = Islet.combined@assays$RNA@counts
sce <- decontX(sce)
Islet.combined@assays$decont = Islet.combined@assays$RNA
Islet.combined@assays$decont@counts = sce$decontXcounts

sce = Islet.combined@assays$RNA@counts
sce <- decontX(sce, background = raw)
#Islet.combined@assays$decont = Islet.combined@assays$RNA
#Islet.combined@assays$decont@data = sce$decontXcounts
Outs =as.data.frame(t(sce$decontXcounts))
DefaultAssay(Islet.combined) = "decont"


umap <- reducedDim(sce, "decontX_UMAP")
plotDimReduceCluster(x = sce$decontX_clusters,
    dim1 = umap[, 1], dim2 = umap[, 2])

DefaultAssay(Beta_Seu) = "RNA"

pdf("Beta_FeaturePlots.pdf", width = 10, height = 10)
print(FeaturePlot(Beta_Seu, unique(c("INS", "DLK1", "SIX3", "PDX1", "CHGA", "ID1", "ID2", "ID3", "MT1F"))))
dev.off()

Set up clustering functions

#GCG     Alpha
#INS     Beta
#SST     Delta
#GHRL    Epsilon
#PPY     Gamma
#PRSS1   Acinar
#KRT19   Ductal
#COL1A1  Stellate
#VWF     Endo
#SDS     Macro
#TPSAB1  Mast
#SOX10   Schwann


GeneLists = list()
GeneLists[["MainList"]] = c("GCG", "INS", "SST", "GHRL", "PPY", "PRSS1", "KRT19", "COL1A1", "VWF", "SDS", "TPSAB1", "SOX10")
GeneLists[["Seger_Alpha1"]] = c("CKS2", "UBE2C", "HMGB2", "PTTG1", "CDKN3", "UBE2T", "CDK1", "CENPW", "TOP2A", "NUSAP1", "PBK", "BIRC5", "CDC20", "CDCA3", "ZWINT", "CDCA8", "AURKB", "KIAA0101", "PRC1", "CCNB1", "TPX2", "HMMR", "ASF1B", "CCNB2", "RRM2", "NDC80", "PLK1", "CENPM", "CCNA2", "BUB1", "CENPF", "TACC3", "TROAP", "FAM64A", "SPAG5", "ANLN", "MK167", "GTSE1", "CENPA", "NEK2", "ASPM", "KIF20A", "CEP55", "DIAPH3", "HJURP", "DEPDC1", "MELK", "KIF23", "CENPE", "DEPDC1B", "IQGAP3CDC25C", "DLGAP5", "CKAP2L", "KIF14", "CIT", "FAM111B")
GeneLists[["Seger_Alpha2"]] = c("PEMT", "C10orf10", "HLA-E", "SLC7A8", "ABCC8", "DSP", "GOPC2", "PPY", "ARRDC4", "PLCE1", "MVP", "ALDOC", "TMEM150A", "SHC2", "GPD1", "C1orf116", "KCNQ10T1", "MEG3", "NPY")
GeneLists[["Seger_Beta"]] = c("ID2", "ID1", "INS", "RBP4", "FFAR4", "ID3")
GeneLists[["Beta_Subpops"]] = c("CFAP126", "DKK3", "PDX1", "GATA2", "CDH1", "NCAM1", "CD9", "ITGA6", "MKI67", "PDGFRA", "MAPK1", "STAT3", "CASP3", "SLC2A2", "DECR1", "CREB1", "AXIN1", "ST8SIA1", "MAFA", "NKX6-1", "PPP1R1A", "ABCC9", "SIX3", "RFX6", "MAFB",  "NEUROD1")
GeneLists[["Seger_Acinar"]] = c("CPB1", "CTRC", "PLA2G1B", "MT1G", "PNLIPRP2", "MT2A", "MT1F", "CELA3A", "MT1X", "MT1E", "FKBP11", "CELAZA", "PNLIPRP1", "TPST2", "SYCN", "REG3G", "PDIA2", "CELA3B", "CLPS", "CTRL", "SLC43A1", "MT1MCEL", "MT1H", "GP2", "AKR7A3", "SLC39A5", "MT1HL1", "TMED6", "SLC38A5", "PRSS2", "SLC30A2", "CXCL12", "CHAC1", "SERPINI2", "CPA4", "FUT1", "GCAT", "LOC100506281", "ENTPD3", "CELA2B", "MT1B", "MAT1A", "SERPINI1", "AMY2B", "FAM46C", "CLPSL1", "MT1A", "MT1L", "AQP8", "PM20D1", "SLC38A3", "IL22RA1", "AQP12B", "KCNE4", "ARHGDIG", "REG4", "DERL3", "AMY2A", "RBPJL", "AQP12A", "COL15A1", "SARDH", "MT1DP", "GPHA2", "CENPM")
GeneLists[["Peng_Ductal"]] = c("FXYD2", "KRT19", "MUC1", "FXYD3", "MKI67", "TOP2A", "CCNB1", "CCNB2", "AMBP")
GeneLists[["Baron_Stellate"]] = c("RGS5", "PDGFRB", "SPARC", "GJA4", "CSPG4", "EDNRB", "COL4A1", "PDGFRA", "COL1A1", "FN1", "THY1", "LUM", "MMP2", "TIMP1", "VCAN", "CXCL3", "IL33", "IL11", "FGF2")
GeneLists[["Rorsman_Pardo_Delta"]] = toupper(c("Crhr2", "Ghsr", "Glp1r", "Gipr", "Chrm3", "Chrm4", "Adra2a", "Sstr1", "Sstr3", "Ffar4", "Insr", "Gcgr", "GABRA1", "GABRB1", "GABRD", "GABRE", "Gabrq", "pdx1", "slc2a2",  "gck"))
GeneLists[["EndoList"]] = c("SLCO1C1", "SLC38A5", "MYH11", "MRC1", "CLDN5", "ITM2A", "FLT1", "VTN",  "COL1A1", "COL3A1", "LUM", "DCN", "PTGDS", "AMBP", "RGS5", "ACTA2", "TAGLN", "ABCC9", "KCNJ8")
GeneLists[["MastCells"]] = c("TPSAB1", "KIT", "MS4A2", "GATA2", "ENPP3", "CD63",  "CD3E", "TRAC", "IL7R", "CCL17", "IDO1", "SDS", "MERTK", "ITGAX")
GeneLists[["GammaCells"]] = c("PPY" , "BTG2", "ID4", "CNTNAP5", "CTD-20008P7.8", "ABCC8", "NEUROD1", "PAX6", "FOLR1", "PYY", "TSPAN8",  "GCG", "INS", "SST", "GHRL","PRSS1", "KRT19" ,"COL1A1" ,"SDS","SKI", "FEV", "FOXA2", "SMAD4", "TCF3", "PHF1", "NR2F6", "FOXO6", "JUND", "ATP1A1", "EPCAM", "LAMP1", "TMEM30B", "DNER", "MGAT4B", "KCNH2", "LRFN4", "NDFIP1", "VEGFA", "PTOV1", "ARX", "LAMP1", "TXNDC5", "CSNK1G2")
GeneLists[["MacroCells"]] = c("CD68", "CD163", "MRC1", "FCGR1", "MAFB", "ITGAX", "FLT3", "LYVE1", "CD3E", "ABCG1", "ITGAM", "ADGRE1", "ICAM2", "TIMD4", "KLF2", "MERTK", "H2-AB1", "CCR2")
GeneLists[["Lori_Extras"]] = c("TTR", "GC", "PCSK2", "GPX3", "SLC7A2", "VGF", "SCG5", "PCSK1NCST3", "IRX1", "IRX2", "FAP", "DPP4", "LOXL4", "INS", "HADH", "PCSK1NKX6-1", "G6PC2", "PAPSS2", "IAPP",  "DLK1", "MAFA", "SIX3", "OLIG1", "ADCYAP1", "PDX1", "GCK", "SST", "POU3F1", "RBP4", "LEPR", "GHRL", "PPY", "ETV1", "RAP1GAP", "PRSS1", "CPA1", "KRT7", "KRT19", "KRT18", "S100A11", "SDC4", "LCN2", "KRT8", "TMSB4X", "ANXA2", "ZFP36L1", "COL1A1", "PDGFRB", "SPARC", "SERPINH1", "TBX3", "CALD1", "COL4A1", "CEBPB", "LGALS1", "EDNRB", "IGFBP7", "PLVAP", "RGCC", "PODXL", "CD93", "ENG", "FLT1", "KDR", "FKBP1A", "SOX18", "VWF", "SDS", "TPSAB1", "SOX10", "CHGA", "CHGB", "SYP", "KCNJ11", "GLUT1", "KRT19", "PRSS1", "HERPUD1", "HSPA5", "DDIT3")





DefaultAssay(Islet.combined) = "RNA"
for(g in names(GeneLists)){
 GeneLists[[g]] = subset(GeneLists[[g]], GeneLists[[g]] %in% row.names(Islet.combined))
}


#Baron et al 2016 - A Single-Cell Transcriptomic Map of the Human and Mouse Pancreas Reveals Inter- and Intra-cell Population Structure - STELLATE MARKERS

DEGs = function(Input){
Input.markers <- FindAllMarkers(Input, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) 
return(Input.markers)}


##Check population on UMAP
#CheckInput = row.names(EmbryonicHypo)
CheckUMAP = function(SeuFilez){
  if(class(CheckInput) == "data.frame"){
CheckMeta = as.data.frame(c(rep("RoI", length(row.names(CheckInput)))))
colnames(CheckMeta) = "Pop"
row.names(CheckMeta) = row.names(CheckInput)
  }else{
    CheckMeta = as.data.frame(c(rep("RoI", length(colnames(CheckInput)))))
colnames(CheckMeta) = "Pop"
row.names(CheckMeta) = colnames(CheckInput)
  }

SeuFilez = AddMetaData(SeuFilez, CheckMeta, "CheckMeta")
Idents(SeuFilez) = "CheckMeta"
DimPlot(SeuFilez, reduction="umap")
}


########################################################################
GenerateMetaData = function(ListMeta){
MetaOutput = as.data.frame(matrix(ncol = 1, nrow =0))  
for(x in seq(1,length(ListMeta),1)){
if(class(ListMeta[[x]]) == "data.frame"){
Temp = as.data.frame(rep(names(ListMeta)[[x]], length(row.names(ListMeta[[x]]))))  
colnames(Temp) = "Pop"
row.names(Temp) = row.names(ListMeta[[x]])
}else{
Temp = as.data.frame(rep(names(ListMeta)[[x]], length(colnames(ListMeta[[x]]))))  
colnames(Temp) = "Pop"
row.names(Temp) = colnames(ListMeta[[x]]) 
}
MetaOutput = rbind(MetaOutput, Temp)  
}
return(MetaOutput)
}
########################################################################



#Fetal.integrated.HypothalV2.Harmony_RNA = as.data.frame(t(as.data.frame(EmbryonicHypo@assays$RNA@counts)))
#All.integrated_DataSubs = Fetal.integrated.HypothalV2.Harmony_RNA %>% dplyr::select(SIM1, `NKX2-1`, MC4R, SST, KISS1, GHRH, TH, TRH, OXT, AVP, CRH, BDNF, PCSK1, GAL,POMC,LEPR, NPY,AGRP, TBX3,  SIX6,  OTP, POU3F2, GABRA5, NPTX2, HCRT, CARTPT, ONECUT2, FOXP2, HDC, CCK, LHX6, LHX9, ARC, NRN1, LHX8, PMCH, `NKX2-2`, LHX1, SIX3, LHX8, PMCH, `NKX2-2`, LHX1, SIX3, SEMA3C, PITX2, PPP1R17, NPY1R, FOXG1, NEUROD6, RGS16,  BARHL1, BCL11A, BCL11B, DMRTA2, EBF2, EMX2, HEY1, LMX1A, MEIS2, MESP1, NHLH2, NHLH1, NME2,  `NKX2-2`, NPAS4, NR2F1, NR2F2, NR5A2, NR5A1, ONECUT1, ONECUT3, OTX2, PBX3, PITX2, PRDM8, SOX14, SOX3, ST18, TBR1, TCF7L2, ZBTB16, ZIC1, ZIC2, LHX1, TAC3, BCL11B, HIVEP2,RAB3B, DNAJC6, DIRAS3, OLFM3, SYT6, S100A10, FAM84A, FAM84A, VSNL1, HS6ST1, CTNNB1, CAMKV, CADM2, KIT, PPP3CA, SPOCK3, DYNC1I1, MAP7D2, PNCK, NEFM, DNAJC12, ATP8A2, KCNMB2)

#All.integrated_DataSubsV2 = Fetal.integrated.HypothalV2.Harmony_RNA %>% select(OligList, AstroMicro, EpendyTanyList, EndoList, MuralPeriVLMCList, BroadGenes, HypothalamicProgenSubclass, WellKnownGenes, SuspectedGenes)


ClusterFunc_All_RNA = function(SeuFile){
Filename = as.character(substitute(SeuFile))

for(y in set.kparam){
for(z in set.dim){
for(v in set.res){
DefaultAssay(SeuFile) = "integrated"
SeuFile <- FindNeighbors(SeuFile, k.param=y, dims=1:z)
SeuFile <- FindClusters(SeuFile, resolution = v)
DefaultAssay(SeuFile) = "RNA"

pdf(paste("ISLETS_", Filename, "_res", v, "_k", y, "_dim", z, "_umapSML.pdf", sep=""), width=12, height=10)
dimplot = DimPlot(SeuFile, reduction="umap", label=T)
print(dimplot)
dev.off()


pdf(paste("ISLETS_", Filename, "_res", v, "_k", y, "_dim", z, "_umapLGE.pdf", sep=""), width=25, height=25)
dimplot = DimPlot(SeuFile, reduction="umap", label=T)
print(dimplot)
dev.off()

###OUTPUT MEAN OF COUNTS
#Subset = subset(All.integrated_DataSubsV2, row.names(All.integrated_DataSubsV2) %in% colnames(SeuFile))
#MetaData = as.data.frame(SeuFile@meta.data$seurat_clusters)
#row.names(MetaData) = row.names(SeuFile@meta.data)
#SubsetwMeta = merge(Subset, MetaData, by=0)
#row.names(SubsetwMeta) = SubsetwMeta$Row.names
#SubsetwMeta$Row.names = NULL

#SubsetSum = Subset %>% group_by(SubsetwMeta$`SeuFile@meta.data$seurat_clusters`) %>% summarise_all(sum, na.rm = TRUE) 
#write.csv(SubsetSum, paste("ISLETS_", Filename, "_SumCounts_k", y, "_dim", z, "_res", v, ".csv", sep=""), row.names = F)

#SubsetMean = Subset %>% group_by(SubsetwMeta$`SeuFile@meta.data$seurat_clusters`) %>% summarise_all(mean, na.rm = TRUE) 
#write.csv(SubsetMean, paste("ISLETS_", Filename, "_MeanCounts_k", y, "_dim", z, "_res", v, ".csv", sep=""), row.names = F)


###VLNS
AllVlns = list()

for(d in names(GeneLists)){
Genes = GeneLists[[d]]
StorePlots = list()
  for(x in Genes[1]){
            plotA <- VlnPlot(SeuFile, features = x, pt.size = 0, same.y.lims = F,)
            plotA <- plotA + coord_flip()+ theme(axis.ticks.x= element_blank(), axis.title.y=element_blank(), axis.text.x=element_blank(), axis.title.x=element_blank(), 
                                                 axis.ticks.y = element_blank(), legend.position = "none", plot.title = element_text(size=12))+ labs(title = d, subtitle = Genes[1])
            StorePlots[[x]] = plotA 
            }
  for(x in Genes[2:length(Genes)]){
            plotB <- VlnPlot(SeuFile, features = x, pt.size = 0, same.y.lims = F,)
            plotB <- plotB + coord_flip()+ theme(axis.ticks.x= element_blank(), axis.title.y=element_blank(), axis.text.x=element_blank(), axis.title.x=element_blank(), 
                      axis.ticks.y = element_blank(), legend.position = "none", axis.text.y = element_blank(), plot.title = element_text(size=12))
           StorePlots[[x]] = plotB
           }
AllVlns[[d]] <- ggarrange(plotlist = StorePlots, widths=c(1.4, rep(1, length(Genes)-1)), ncol = 60,  nrow = 1)  
}
pdf(paste("ISLETS_", Filename, "_res", v, "_k", y, "_dim", z, "_AllMultiVlns.pdf", sep=""), width=60, height=length(unique(SeuFile@active.ident)))
print(AllVlns)
dev.off()

#Cell No
CellNo = as.data.frame(table(SeuFile@meta.data$seurat_clusters))
write.csv(CellNo, paste("ISLETS_", Filename, "_counts_k", y, "_dim", z, "_res", v, ".csv", sep=""), row.names = F)
}}

#Feature Plots  
FPList = list()  

for(d in names(GeneLists)){
Genes = GeneLists[[d]]

FPSinglePage = list()
FPSinglePage[[1]] = FeaturePlot(SeuFile, Genes[1], reduction="umap") + labs(title=paste(d, Genes[1]))
for(p in seq(2, length(Genes), 1)){
FPSinglePage[[p]] = FeaturePlot(SeuFile, Genes[p], reduction="umap") 
}
FPList[[d]] = ggarrange(plotlist = FPSinglePage, ncol=6, nrow=10)
}

pdf(paste("ISLETS_", Filename, "_dim", z, "_AllFPs.pdf", sep=""), width = 25, height = 45)
print(FPList)
dev.off()
}}

Quality control

Islets.Data= list()
QCViolins = list()
QCHist = list()
QCTable = as.data.frame(matrix(ncol=10, nrow=0))


for(x in c("HP17225", "HP19345", "HP20217", "HP17117", "HP18320", "HP19208", "HP18227", "HP18047", "HP17209", "HP17250", "HP17082")){
Islets.Data[[x]] = Read10X(data.dir = paste("~/Dropbox/Columbia/Collins Lab Collab/Collins Islet Data/islet_cellranger_filtered/", x, "_filtered_feature_bc_matrix", sep=""))
ScreeingFile = CreateSeuratObject(counts = Islets.Data[[x]], project = x, min.cells = 3, min.features = 0)
ScreeingFile[["percent.mt"]] <- PercentageFeatureSet(ScreeingFile, pattern = "^MT-")
QCViolins[[x]] = VlnPlot(ScreeingFile, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
TableData = as.data.frame(t(c(x, median(ScreeingFile$nFeature_RNA), mean(ScreeingFile$nFeature_RNA), sd(ScreeingFile$nFeature_RNA), median(ScreeingFile$nCount_RNA), mean(ScreeingFile$nCount_RNA), sd(ScreeingFile$nCount_RNA), median(ScreeingFile$percent.mt), mean(ScreeingFile$percent.mt), sd(ScreeingFile$percent.mt))))
colnames(TableData) = c("Sample", "nFeature_Median", "nFeature_Mean", "nFeature_SD", "nCount_Median", "nCount_Mean", "nCount_SD", "percent.mt_Median", "percent.mt_Mean", "percent.mt_SD")
QCTable = rbind(QCTable, TableData)
QCHist[[x]] = ggplot(ScreeingFile@meta.data, aes(x=nFeature_RNA)) + 
  geom_histogram(binwidth=30)+ ggtitle(x)+ scale_x_continuous(breaks = seq(0,max(ScreeingFile@meta.data$nFeature_RNA),200)) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 
}


pdf(paste("ISLETS_QC_Violins_24MAR22.pdf", sep=""), width=20, height=8)
print(QCViolins)
dev.off()

pdf(paste("ISLETS_QC_Histograms_24MAR22.pdf", sep=""), width=15, height=8)
print(QCHist)
dev.off()

Integration and main level clustering

Islets.DataV3 = Islets.Data
Islets.DataV3[["HP17117"]] = NULL
Islets.DataV3[["HP18320"]] = NULL
Islets.DataV3[["HP17250"]] = NULL
Islets.DataV3[["HP17209"]] = NULL
Islets.DataV3[["HP17225"]] = NULL
Islets.DataV4 = list()

for(x in 1:length(Islets.DataV3)){
File =CreateSeuratObject(counts = Islets.DataV3[[x]], project = names(Islets.DataV3)[x], min.cells = 3, min.features = 1000)
Temp = as.data.frame(File@assays$RNA@counts)
colnames(Temp) = paste0(names(Islets.DataV3)[x], colnames(Temp), sep="_")
Temp =CreateSeuratObject(counts = Temp, project = names(Islets.DataV3)[x], min.cells = 3, min.features = 0)
Temp <- NormalizeData(Temp)
Temp <- FindVariableFeatures(Temp, do.plot = F, display.progress = F)

Temp = ScaleData(Temp, vars.to.regress = c("nFeature_RNA", "percent_mito"),
    verbose = F)
Temp = RunPCA(Temp, verbose = F, npcs = 20)
Temp = RunUMAP(Temp, dims = 1:10, verbose = F)

nExp_poi <- round(0.1*nrow(Temp@meta.data))
Temp<- doubletFinder_v3(Temp, pN = 0.25, pK = 0.09, nExp = nExp_poi, PCs = 1:10)
Temp.meta = Temp@meta.data
Temp.meta2 = subset(Temp.meta, Temp.meta[[length(Temp.meta)]] == "Singlet")
Islets.DataV4[[x]] = subset(Temp, cells = row.names(Temp.meta2))
}

features <- SelectIntegrationFeatures(object.list = Islets.DataV4)
Islet.anchors <- FindIntegrationAnchors(object.list = Islets.DataV4, anchor.features = features)
Islet.combined <- IntegrateData(anchorset = Islet.anchors)
DefaultAssay(Islet.combined) <- "integrated"

save(list=c("Islet.combined"), file = "~/Islet_StdIntegration_DbltRm.RData")
DefaultAssay(Islet.combined) <- "RNA"

MT.genes <- grep(pattern = "^MT-", x = rownames(x = Islet.combined), value = TRUE)
DefaultAssay(Islet.combined) <- "integrated"

percent.mt <- Matrix::colSums(Islet.combined@assays[["RNA"]][MT.genes, ])/Matrix::colSums(Islet.combined@assays[["RNA"]])
Islet.combined  = AddMetaData(Islet.combined, percent.mt, "percent.mt")
Islet.combined <- ScaleData(Islet.combined, vars.to.regress = "percent.mt")
Islet.combined <- RunPCA(Islet.combined, npcs = 30, verbose = FALSE)
Islet.combined <- RunUMAP(Islet.combined, reduction = "pca", dims = 1:30)
save(list=c("Islet.combined"), file = "~/Islet_StdIntegration_DbltRm.RData")

DPlist2 = list()
DPlist = list()
FPlist = list()
Filename = "Islets_DbltsRm_1JUL22"
for(y in seq(10,50,10)){
DefaultAssay(Islet.combined) = "integrated"
Islet.combined = RunPCA(Islet.combined, npcs = y)
for(z in seq(10,y,10)){
for(s in c(2,5,10)){
Islet.combined <- RunUMAP(Islet.combined, dims = 1:z, spread= s)
DefaultAssay(Islet.combined) = "RNA"

FPlist[[paste("PCA", y, "_dims", z, "_spread", s)]] = FeaturePlot(Islet.combined, c("GCG", "INS", "SST", "GHRL", "PPY", "PRSS1", "KRT19", "COL1A1", "VWF", "SDS", "TPSAB1", "SOX10"), reduction="umap") 

DPlist[[paste("PCA", y, "_dims", z, "_spread", s)]] = DimPlot(Islet.combined, reduction="umap", split.by = "orig.ident") + labs(title = paste("PCA", y, "_dims", z, "_spread", s))

DPlist2[[paste("PCA", y, "_dims", z, "_spread", s)]] = DimPlot(Islet.combined, reduction="umap", group.by = "orig.ident", label=T) + labs(title = paste("PCA", y, "dims", z, "spread", s))
}}}  
pdf(paste("ISLETS_", Filename, "_FeaturePlots.pdf", sep=""), width=20, height=20)
print(FPlist)
dev.off()

pdf(paste("ISLETS_", Filename, "_SPLITUMAP.pdf", sep=""), width=50, height=5)
print(DPlist)
dev.off()

pdf(paste("ISLETS_", Filename, "_GROUPUMAP.pdf", sep=""), width=8, height=5)
print(DPlist2)
dev.off()



DefaultAssay(Islet.combined) = "integrated"
Islet.combined <- RunPCA(Islet.combined, npcs = 20, verbose = FALSE)
Islet.combined <- RunUMAP(Islet.combined, reduction = "pca", dims = 1:20, spread=2)
Islet.combined_UMAP = as.data.frame(Islet.combined@reductions$umap@cell.embeddings)

Alpha_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > -5 &  Islet.combined_UMAP$UMAP_1 < 10 & Islet.combined_UMAP$UMAP_2 > 5)
Epsilon_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > 13 &  Islet.combined_UMAP$UMAP_1 < 15.7 & Islet.combined_UMAP$UMAP_2 > 2 & Islet.combined_UMAP$UMAP_2 < 3)
Gamma_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > 12  & Islet.combined_UMAP$UMAP_2 > 2 & ! row.names(Islet.combined_UMAP) %in% row.names(Epsilon_Clust))
Delta_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > 5 &  Islet.combined_UMAP$UMAP_1 < 14 & Islet.combined_UMAP$UMAP_2 > -3 & Islet.combined_UMAP$UMAP_2 < 3)
Acinar_Clust = subset(Islet.combined_UMAP,  Islet.combined_UMAP$UMAP_1 < -9 & Islet.combined_UMAP$UMAP_2 > -2)
Beta_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > -3 &  Islet.combined_UMAP$UMAP_1 < 10 & Islet.combined_UMAP$UMAP_2 < -3)
Cycling_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > 12  & Islet.combined_UMAP$UMAP_2 < -3)
Stellate_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > -12 &  Islet.combined_UMAP$UMAP_1 < -5 & Islet.combined_UMAP$UMAP_2 < -11.6 & Islet.combined_UMAP$UMAP_2 > -15)
Endo_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 > -12 &  Islet.combined_UMAP$UMAP_1 < -5 & Islet.combined_UMAP$UMAP_2 < -9.5 & Islet.combined_UMAP$UMAP_2 > -11.6)
Ductal_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 < -5 & Islet.combined_UMAP$UMAP_2 < -2  & Islet.combined_UMAP$UMAP_2 > -9.5)
Macro_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 < -5  & Islet.combined_UMAP$UMAP_2 < -17.7)
WBC_Clust = subset(Islet.combined_UMAP, Islet.combined_UMAP$UMAP_1 < -5 & Islet.combined_UMAP$UMAP_2 < -15 & Islet.combined_UMAP$UMAP_2 > -17.7)

CheckInput = Acinar_Clust
CheckUMAP(Islet.combined)

#Generate Metadata
Islets_Barcs_List = list("Beta" = Beta_Clust, "Acinar" = Acinar_Clust, "Ductal" = Ductal_Clust, "Alpha" = Alpha_Clust, "Macro" = Macro_Clust,  "Delta" = Delta_Clust, "Gamma" = Gamma_Clust, "Epsilon" = Epsilon_Clust, "Cycling" = Cycling_Clust, "WBCs" = WBC_Clust, "Endo" = Endo_Clust,"Stellate" = Stellate_Clust)
Islets_Barcs = GenerateMetaData(Islets_Barcs_List)
Islet.combined = AddMetaData(Islet.combined, Islets_Barcs, "Islets_Barcs")
Idents(Islet.combined) = "Islets_Barcs"

pdf("Islets_DbltsRm_GHRL.pdf", width = 16, height=16)
print(FeaturePlot(Islet.combined, "GHRL"))
dev.off()

pdf("Islets_DbltsRm_UMAP.pdf", width = 16, height=16)
print(DimPlot(Islet.combined, label = T) + NoLegend())
dev.off()
DefaultAssay(Islet.combined) = "RNA"
pdf("Islets_DbltsRm_FeaturePlots.pdf", width = 15, height=10)
print(FeaturePlot(Islet.combined, c("GCG", "INS", "SST", "GHRL", "PPY", "PRSS1", "KRT19", "COL1A1", "VWF", "SDS", "TPSAB1", "SOX10"), reduction="umap"))
dev.off()

Check overlap with Azimuth

### Output data to check Azimuth labels
saveRDS(Islet.combined, file = "Islet.combinedDbltRm.rds")

predictions <- read.delim('~/Downloads/azimuth_pred_DbltRm.tsv', row.names = 1)
Islets_Barcs2 = Islets_Barcs
Islets_Barcs2$Pop = gsub("WBCs", "Immune", gsub("Macro", "Immune", Islets_Barcs2$Pop))
predictions$predicted.annotation.l1.score = NULL
predictions$mapping.score = NULL
predictions$predicted.annotation.l1 = gsub("beta", "Beta", gsub("delta", "Delta", gsub("ductal", "Ductal", gsub("alpha", "Alpha", gsub("gamma", "Gamma", gsub("immune", "Immune", gsub("epsilon", "Epsilon", gsub("acinar", "Acinar", gsub("activated_stellate", "Stellate", gsub("endothelial", "Endo", gsub("quiescent_stellate", "Stellate", gsub("schwann", "Schwann", gsub("cycling", "Cycling", predictions$predicted.annotation.l1)))))))))))))


Merge_Barcs = merge(Islets_Barcs2, predictions, by = 0)
Merge_Barcs$Equal = Merge_Barcs$Pop == Merge_Barcs$predicted.annotation.l1

Merge_False = subset(Merge_Barcs, Merge_Barcs$Equal == F)
Merge_False$Overlap = paste(Merge_False$Pop, Merge_False$predicted.annotation.l1, sep="_")
OverlapTable = as.data.frame(table(Merge_False$Overlap))
AB = subset(Merge_False, Merge_False$Overlap == c("Beta_Alpha", "Beta_Delta"))
row.names(AB) = AB$Row.names
AB$Row.names = NULL
AB$Pop = NULL
AB$predicted.annotation.l1 = NULL
AB$Equal = NULL
Islet.combined = AddMetaData(Islet.combined, AB, "Beta_Alpha")
Idents(Islet.combined) = "Beta_Alpha"
DP = DimPlot(Islet.combined, label = T, repel=T)+ NoLegend() 

pdf("Islets_Discreprant.pdf", width = 8, height = 8)
print(DP)
dev.off()

Subclustering data

### Iterate dimensions for each of the subclusters
Idents(Islet.combined) = "Islets_Barcs"
for(b in c("Stellate", "Endo")){ #unique(Islets_Barcs$Pop)
SubsetSeu = subset(Islet.combined, idents = b)

if(b == "Stellate"){PlotGenes = c("COL1A1", "COL1A1", "PDGFRB", "SPARC", "SERPINH1", "TBX3", "CALD1", "COL4A1", "CEBPB", "LGALS1", "EDNRB", "IGFBP7", GeneLists[["Baron_Stellate"]])}else if(b == "Endo"){PlotGenes = c("CSPG4", "GJA4", "APLN", "CD82", "CHST1", "PTH1R", "BMX", "GKN3", "IGFBP2", "DCN", "LIM", "PDGFRA", "IL33", "PLVAP", "RGCC", "PODXL", "CD93", "ENG", "FLT1", "KDR", "FKBP1A", "SOX18", "VWF", GeneLists[["EndoList"]])}else if(b == "Beta"){PlotGenes = c("INS", "INS", "HADH", "PCSK1", "NKX6-1", "G6PC2", "PAPSS2", "IAPP",  "DLK1", "MAFA", "SIX3", "OLIG1", "ADCYAP1", "PDX1", "GCK", "CHGA", "CHGB", "SYP", "KCNJ11", GeneLists[["Seger_Beta"]], GeneLists[["Beta_Subpops"]])}else if(b == "Acinar"){PlotGenes = c("PRSS1","PRSS1", "CPA1", GeneLists[["Seger_Acinar"]])}else if(b == "Ductal"){PlotGenes = c("KRT7", "KRT19", "KRT18", "S100A11", "SDC4", "LCN2", "KRT8", "TMSB4X", "ANXA2", "ZFP36L1", GeneLists[["Peng_Ductal"]])}else if(b == "Alpha"){PlotGenes = c("GCG", "TTR", "GC", "PCSK2", "GPX3", "SLC7A2", "VGF", "SCG5", "PCSK1NCST3", "IRX1", "IRX2", "FAP", "DPP4", "LOXL4", GeneLists[["Seger_Alpha1"]], GeneLists[["Seger_Alpha2"]])}else if(b == "Macro"){PlotGenes = c(GeneLists[["MacroCells"]])}else if(b == "WBCs"){PlotGenes = c(GeneLists[["MainList"]], "TPSAB1", GeneLists[["MastCells"]])}else if(b == "Delta"){PlotGenes = c("SST","SST", "POU3F1", "RBP4", "LEPR", GeneLists[["Rorsman_Pardo_Delta"]])}else if(b == "Gamma"){PlotGenes = c("PPY" ,"PPY", "ETV1", "RAP1GAP",  GeneLists[["GammaCells"]]
)}else{PlotGenes = c(GeneLists[["MainList"]])}
PlotGenes = unique(PlotGenes)

DPlist2 = list()
DPlist = list()
FPlist = list()
Filename = paste(b, "27JUL22", sep="")
for(y in seq(10,50,10)){
DefaultAssay(SubsetSeu) = "integrated"
SubsetSeu = RunPCA(SubsetSeu, npcs = y)
for(z in seq(10,y,10)){
for(s in c(2,5,10)){
SubsetSeu <- RunUMAP(SubsetSeu, dims = 1:z, spread= s)
DefaultAssay(SubsetSeu) = "RNA"

FPlist[[paste("PCA", y, "_dims", z, "_spread", s)]] = FeaturePlot(SubsetSeu, PlotGenes, reduction="umap") 

DPlist[[paste("PCA", y, "_dims", z, "_spread", s)]] = DimPlot(SubsetSeu, reduction="umap", split.by = "orig.ident") + labs(title = paste("PCA", y, "_dims", z, "_spread", s))

DPlist2[[paste("PCA", y, "_dims", z, "_spread", s)]] = DimPlot(SubsetSeu, reduction="umap", group.by = "orig.ident", label=T) + labs(title = paste("PCA", y, "dims", z, "spread", s))
}}}  
pdf(paste("ISLETS_", Filename, "_FeaturePlots.pdf", sep=""), width=20, height=length(PlotGenes)/1.3)
print(FPlist)
dev.off()

pdf(paste("ISLETS_", Filename, "_SPLITUMAP.pdf", sep=""), width=50, height=5)
print(DPlist)
dev.off()

pdf(paste("ISLETS_", Filename, "_GROUPUMAP.pdf", sep=""), width=8, height=5)
print(DPlist2)
dev.off()
}


### SUBCLUSTER ALPHA
Idents(Islet.combined) = "Islets_Barcs"
Alpha_Seu = subset(Islet.combined, idents = "Alpha")
DefaultAssay(Alpha_Seu) = "integrated"
Alpha_Seu <- RunPCA(Alpha_Seu, npcs = 20, verbose = FALSE)
Alpha_Seu <- RunUMAP(Alpha_Seu, reduction = "pca", dims = 1:20, spread=2)
Alpha_Seu_UMAP = as.data.frame(Alpha_Seu@reductions$umap@cell.embeddings)

set.res = c(1)
set.dim = c(10)
set.kparam = c(10)
#ClusterFunc_All_RNA(Alpha_Seu)

DefaultAssay(Alpha_Seu) = "integrated"
Alpha_Seu <- FindNeighbors(Alpha_Seu, k.param=10, dims=1:10)
Alpha_Seu <- FindClusters(Alpha_Seu, resolution = 1)
Alpha_Clust1_ToClean = subset(Alpha_Seu, idents = c(5,13))
Alpha_Clust2_ToClean = subset(Alpha_Seu, idents = c(10))

Alpha_Clust1 = subset(Alpha_Seu_UMAP, row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Clust1_ToClean) &  Alpha_Seu_UMAP$UMAP_1 > 0 & Alpha_Seu_UMAP$UMAP_2 < -3  | Alpha_Seu_UMAP$UMAP_1 > -1 & Alpha_Seu_UMAP$UMAP_2 < -5 )

Alpha_Clust2 = subset(Alpha_Seu_UMAP, row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Clust2_ToClean) &  Alpha_Seu_UMAP$UMAP_1 < 0 & Alpha_Seu_UMAP$UMAP_2 < -2  | Alpha_Seu_UMAP$UMAP_1 < -1 & Alpha_Seu_UMAP$UMAP_2 < -4 )

Alpha_Rest = subset(Alpha_Seu, cells = c(row.names(Alpha_Clust1), row.names(Alpha_Clust2)), invert=T)

set.res = c(1)
set.dim = c(10)
set.kparam = c(5,4,3)
#ClusterFunc_All_RNA(Alpha_Rest)

DefaultAssay(Alpha_Rest) = "integrated"
Alpha_Rest <- FindNeighbors(Alpha_Rest, k.param=5, dims=1:10)
Alpha_Rest <- FindClusters(Alpha_Rest, resolution = 1)
Alpha_Clust3_ToClean = subset(Alpha_Rest, idents = c(2,3,6,14,9,10))

Alpha_Clust3 = subset(Alpha_Seu_UMAP, row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Clust3_ToClean) &  Alpha_Seu_UMAP$UMAP_1 > -3 | row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Rest) &  Alpha_Seu_UMAP$UMAP_1 > 3 )
CheckInput = Alpha_Clust3
CheckUMAP(Alpha_Rest)

Alpha_Rest2 = subset(Alpha_Seu, cells = c(row.names(Alpha_Clust1), row.names(Alpha_Clust2), row.names(Alpha_Clust3)), invert=T)

set.res = c(1)
set.dim = c(20)
set.kparam = c(20,50)
ClusterFunc_All_RNA(Alpha_Rest2)

DefaultAssay(Alpha_Rest2) = "integrated"
Alpha_Rest2 <- FindNeighbors(Alpha_Rest2, k.param=20, dims=1:20)
Alpha_Rest2 <- FindClusters(Alpha_Rest2, resolution = 1)
Alpha_Clust4_ToClean = subset(Alpha_Rest2, idents = c(0,8,5))
Alpha_Clust5_ToClean = subset(Alpha_Rest2, idents = c(6))

CheckInput =Alpha_Clust5_ToClean
CheckUMAP(Alpha_Rest2)

Alpha_Clust5 = subset(Alpha_Seu_UMAP, row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Clust5_ToClean) &  Alpha_Seu_UMAP$UMAP_1 < -1 & Alpha_Seu_UMAP$UMAP_1 > -5 | row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Rest2) &  Alpha_Seu_UMAP$UMAP_1 > -4 & Alpha_Seu_UMAP$UMAP_1 < -2.5 & Alpha_Seu_UMAP$UMAP_2 > -1 & Alpha_Seu_UMAP$UMAP_2 < 1)

Alpha_Clust4 = subset(Alpha_Seu_UMAP, row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Clust4_ToClean) &  Alpha_Seu_UMAP$UMAP_1 > -4.5 & Alpha_Seu_UMAP$UMAP_2 < 4 & ! row.names(Alpha_Seu_UMAP) %in% row.names(Alpha_Clust5) | row.names(Alpha_Seu_UMAP) %in% colnames(Alpha_Rest2) &  Alpha_Seu_UMAP$UMAP_1 > -1 & Alpha_Seu_UMAP$UMAP_2 < 1 & ! row.names(Alpha_Seu_UMAP) %in% row.names(Alpha_Clust5))


Alpha_Clust6 = subset(Alpha_Seu, cells = c(row.names(Alpha_Clust1), row.names(Alpha_Clust2), row.names(Alpha_Clust3), row.names(Alpha_Clust4), row.names(Alpha_Clust5)), invert=T)

CheckInput =Alpha_Clust6
CheckUMAP(Alpha_Rest2)

Alpha_Barcs_List = list("Alpha_1" = Alpha_Clust1, "Alpha_2" = Alpha_Clust2, "Alpha_3" = Alpha_Clust3, "Alpha_4" = Alpha_Clust4, "Alpha_5" = Alpha_Clust5, "Alpha_6" = Alpha_Clust6)
Alpha_Barcs = GenerateMetaData(Alpha_Barcs_List)
###



### SUBCLUSTER BETA #Use 16
Beta_Seu = subset(Islet.combined, idents = "Beta")
DefaultAssay(Beta_Seu) = "integrated"
Beta_Seu <- RunPCA(Beta_Seu, npcs = 30, verbose = FALSE)
Beta_Seu <- RunUMAP(Beta_Seu, reduction = "pca", dims = 1:30, spread=2)
Beta_Seu_UMAP = as.data.frame(Beta_Seu@reductions$umap@cell.embeddings)

Beta_Clust1 = subset(Beta_Seu_UMAP,  Beta_Seu_UMAP$UMAP_1 > 1.5 & Beta_Seu_UMAP$UMAP_2 < -4 | Beta_Seu_UMAP$UMAP_1 > 2.5 & Beta_Seu_UMAP$UMAP_2 < -3)

CheckInput = Beta_Clust1
CheckUMAP(Beta_Seu)

Beta_Seu_Rest = subset(Beta_Seu, cells = row.names(Beta_Clust1), invert=T)

set.res = c(1)
set.dim = c(30)
set.kparam = c(50)
ClusterFunc_All_RNA(Beta_Seu)

DefaultAssay(Beta_Seu_Rest) = "integrated"
Beta_Seu_Rest <- FindNeighbors(Beta_Seu_Rest, k.param=50, dims=1:30)
Beta_Seu_Rest <- FindClusters(Beta_Seu_Rest, resolution = 1)
Beta_Clust2 = subset(Beta_Seu_Rest, idents = c(5))
Beta_Clust3 = subset(Beta_Seu_Rest, idents = c(6))

CheckInput = Beta_Clust3_Clean
CheckUMAP(Beta_Seu_Rest)

Beta_Clust2_Clean = subset(Beta_Seu_UMAP, row.names(Beta_Seu_UMAP) %in% colnames(Beta_Clust2) & Beta_Seu_UMAP$UMAP_1 < -5 & Beta_Seu_UMAP$UMAP_2 < 3)

Beta_Clust3_Clean = subset(Beta_Seu_UMAP, row.names(Beta_Seu_UMAP) %in% colnames(Beta_Clust3) & Beta_Seu_UMAP$UMAP_1 < 0 & Beta_Seu_UMAP$UMAP_2 < -2 | row.names(Beta_Seu_UMAP) %in% colnames(Beta_Seu_Rest) & Beta_Seu_UMAP$UMAP_1 < -3 & Beta_Seu_UMAP$UMAP_2 < -4)

Beta_Seu_Rest2 = subset(Beta_Seu, cells = c(row.names(Beta_Clust1), row.names(Beta_Clust2_Clean), row.names(Beta_Clust3_Clean)), invert=T)

set.res = c(1)
set.dim = c(30)
set.kparam = c(5)
#ClusterFunc_All_RNA(Beta_Seu_Rest2)

DefaultAssay(Beta_Seu_Rest2) = "integrated"
Beta_Seu_Rest2 <- FindNeighbors(Beta_Seu_Rest2, k.param=200, dims=1:30)
Beta_Seu_Rest2 <- FindClusters(Beta_Seu_Rest2, resolution = 1)
Beta_Clust4 = subset(Beta_Seu_Rest2, idents = c(0))

Beta_Clust4_Clean = subset(Beta_Seu_UMAP, row.names(Beta_Seu_UMAP) %in% colnames(Beta_Clust4) & Beta_Seu_UMAP$UMAP_1 < 1 & Beta_Seu_UMAP$UMAP_2 > -2 | row.names(Beta_Seu_UMAP) %in% colnames(Beta_Seu_Rest2) & Beta_Seu_UMAP$UMAP_1 < -2 & Beta_Seu_UMAP$UMAP_2 > 0)

Beta_Seu_Rest3 = subset(Beta_Seu, cells = c(row.names(Beta_Clust1), row.names(Beta_Clust2_Clean), row.names(Beta_Clust3_Clean), row.names(Beta_Clust4_Clean)), invert=T)

DefaultAssay(Beta_Seu_Rest3) = "integrated"
Beta_Seu_Rest3 <- FindNeighbors(Beta_Seu_Rest3, k.param=3, dims=1:10)
Beta_Seu_Rest3 <- FindClusters(Beta_Seu_Rest3, resolution = 1)
Beta_Clust4Pt2 = subset(Beta_Seu_Rest3, idents = c(11,13,14,15,28))

Beta_Seu_Rest4 = subset(Beta_Seu, cells = c(row.names(Beta_Clust1), row.names(Beta_Clust2_Clean), row.names(Beta_Clust3_Clean), row.names(Beta_Clust4_Clean), colnames(Beta_Clust4Pt2)), invert=T)

set.res = c(1)
set.dim = c(10,20)
set.kparam = c(5)
#ClusterFunc_All_RNA(Beta_Seu_Rest4)

DefaultAssay(Beta_Seu_Rest4) = "integrated"
Beta_Seu_Rest4 <- FindNeighbors(Beta_Seu_Rest4, k.param=5, dims=1:20)
Beta_Seu_Rest4 <- FindClusters(Beta_Seu_Rest4, resolution = 1)
Beta_Clust5 = subset(Beta_Seu_Rest4, idents = c(3,7))
Beta_Clust6 = subset(Beta_Seu_Rest4, idents = c(1,9))
Beta_Clust7 = subset(Beta_Seu_Rest4, idents = c(8))

set.res = c(1)
set.dim = c(30)
set.kparam = c(10)
#ClusterFunc_All_RNA(Beta_Clust5)

Beta_Clust5 <- FindNeighbors(Beta_Clust5, k.param=10, dims=1:30)
Beta_Clust5 <- FindClusters(Beta_Clust5, resolution = 1)
Beta_Clust5_V2 = subset(Beta_Clust5, idents = c(1))
Beta_Clust5_Clean = subset(Beta_Seu_UMAP, row.names(Beta_Seu_UMAP) %in% colnames(Beta_Clust5_V2) & Beta_Seu_UMAP$UMAP_2 < -2.5 | row.names(Beta_Seu_UMAP) %in% colnames(Beta_Seu_Rest4) & Beta_Seu_UMAP$UMAP_2 < -4)
#Beta_Clust8= subset(Beta_Clust5, cells = row.names(Beta_Clust5_Clean), invert=T)

Beta_Clust6_Clean = subset(Beta_Seu_UMAP, row.names(Beta_Seu_UMAP) %in% colnames(Beta_Clust6) & Beta_Seu_UMAP$UMAP_2 < 4 & Beta_Seu_UMAP$UMAP_1 > 2 | row.names(Beta_Seu_UMAP) %in% colnames(Beta_Seu_Rest4) & Beta_Seu_UMAP$UMAP_1 > 7)

Beta_Clust7_Clean = subset(Beta_Seu_UMAP, row.names(Beta_Seu_UMAP) %in% colnames(Beta_Clust7) | row.names(Beta_Seu_UMAP) %in% colnames(Beta_Seu_Rest4) & Beta_Seu_UMAP$UMAP_2 > 4.5)


Beta_Seu_Rest5 = subset(Beta_Seu, cells = c(row.names(Beta_Clust1), row.names(Beta_Clust2_Clean), row.names(Beta_Clust3_Clean), row.names(Beta_Clust4_Clean), colnames(Beta_Clust4Pt2), row.names(Beta_Clust5_Clean), row.names(Beta_Clust6_Clean), row.names(Beta_Clust7_Clean)), invert=T)


CheckInput = Beta_Clust8
CheckUMAP(Beta_Seu_Rest5)

set.res = c(1)
set.dim = c(20,10)
set.kparam = c(50)
#ClusterFunc_All_RNA(Beta_Seu_Rest5)

Beta_Seu_Rest5 <- FindNeighbors(Beta_Seu_Rest5, k.param=50, dims=1:30)
Beta_Seu_Rest5 <- FindClusters(Beta_Seu_Rest5, resolution = 1)
Beta_Clust8 = subset(Beta_Seu_Rest5, idents = c(2))

Beta_Clust8_Clean = subset(Beta_Seu_UMAP, row.names(Beta_Seu_UMAP) %in% colnames(Beta_Clust8) & Beta_Seu_UMAP$UMAP_2 < 1 & Beta_Seu_UMAP$UMAP_1 > -2 )
Beta_Clust9 = subset(Beta_Seu_Rest5, cells = row.names(Beta_Clust8_Clean), invert=T)

CheckInput = Beta_Clust8_Clean
CheckUMAP(Beta_Seu_Rest5)


Beta_Barcs_List = list("Beta_1" = Beta_Clust1, "Beta_2" = Beta_Clust2_Clean, "Beta_3" = Beta_Clust3_Clean, "Beta_4" = Beta_Clust4_Clean, "Beta_4" = Beta_Clust4Pt2, "Beta_5" = Beta_Clust5_Clean, "Beta_6" = Beta_Clust6_Clean, "Beta_7" = Beta_Clust7_Clean, "Beta_8" = Beta_Clust8_Clean, "Beta_9" = Beta_Clust9)
Beta_Barcs = GenerateMetaData(Beta_Barcs_List)
###


### SUBCLUSTER ACINAR # Use 3
Idents(Islet.combined) = "Islets_Barcs"
Acinar_Seu = subset(Islet.combined, idents = "Acinar")
DefaultAssay(Acinar_Seu) = "integrated"
Acinar_Seu <- RunPCA(Acinar_Seu, npcs = 10, verbose = FALSE)
Acinar_Seu <- RunUMAP(Acinar_Seu, reduction = "pca", dims = 1:10, spread=10)
Acinar_Seu_UMAP = as.data.frame(Acinar_Seu@reductions$umap@cell.embeddings)

set.res = c(1)
set.dim = c(10)
set.kparam = c(50)
#ClusterFunc_All_RNA(Acinar_Seu)

DefaultAssay(Acinar_Seu) = "integrated"
Acinar_Seu <- FindNeighbors(Acinar_Seu, k.param=50, dims=1:10)
Acinar_Seu <- FindClusters(Acinar_Seu, resolution = 1)
Acinar_Clust1 = subset(Acinar_Seu, idents = c(6))

Acinar_Clust1_Clean = subset(Acinar_Seu_UMAP, row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Clust1) & Acinar_Seu_UMAP$UMAP_2 < -5 | Acinar_Seu_UMAP$UMAP_2 < -20 & Acinar_Seu_UMAP$UMAP_1 > 15)

CheckInput = Acinar_Clust1_Clean
CheckUMAP(Acinar_Seu)

Acinar_Rest = subset(Acinar_Seu, cells = row.names(Acinar_Clust1_Clean), invert=T)

set.res = c(1)
set.dim = c(10)
set.kparam = c(5)
#ClusterFunc_All_RNA(Acinar_Rest)

DefaultAssay(Acinar_Rest) = "integrated"
Acinar_Rest <- FindNeighbors(Acinar_Rest, k.param=5, dims=1:10)
Acinar_Rest <- FindClusters(Acinar_Rest, resolution = 1)
Acinar_Clust2_ToClean = subset(Acinar_Rest, idents = c(0))
Acinar_Clust3_ToClean = subset(Acinar_Rest, idents = c(2,6,11))


Acinar_Clust2 = subset(Acinar_Seu_UMAP, row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Clust2_ToClean) & Acinar_Seu_UMAP$UMAP_1> 2 & Acinar_Seu_UMAP$UMAP_2 < 10 | row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Rest) & Acinar_Seu_UMAP$UMAP_1> 15 & Acinar_Seu_UMAP$UMAP_2 < 2 )

Acinar_Clust3 = subset(Acinar_Seu_UMAP, row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Clust3_ToClean) & Acinar_Seu_UMAP$UMAP_1 > 12 & Acinar_Seu_UMAP$UMAP_2 > -3  & ! row.names(Acinar_Seu_UMAP) %in% row.names(Acinar_Clust2) | Acinar_Seu_UMAP$UMAP_1 > 18 & Acinar_Seu_UMAP$UMAP_2 > -2 & ! row.names(Acinar_Seu_UMAP) %in% c(row.names(Acinar_Clust2), row.names(Acinar_Clust1_Clean)))


CheckInput = Acinar_Clust3
CheckUMAP(Acinar_Seu)



Acinar_Rest2 = subset(Acinar_Seu, cells = c(row.names(Acinar_Clust1_Clean), row.names(Acinar_Clust2), row.names(Acinar_Clust3)), invert=T)

set.res = c(1)
set.dim = c(10)
set.kparam = c(30)
ClusterFunc_All_RNA(Acinar_Rest2)

Acinar_Rest2 <- FindNeighbors(Acinar_Rest2, k.param=30, dims=1:10)
Acinar_Rest2 <- FindClusters(Acinar_Rest2, resolution = 1)
Acinar_Clust4_ToClean = subset(Acinar_Rest2, idents = c(2))
Acinar_Clust5_ToClean = subset(Acinar_Rest2, idents = c(0))
Acinar_Clust6_ToClean = subset(Acinar_Rest2, idents = c(3))
Acinar_Clust7_ToClean = subset(Acinar_Rest2, idents = c(4,5))
Acinar_Clust8_ToClean = subset(Acinar_Rest2, idents = c(1))

Acinar_Clust4 = subset(Acinar_Seu_UMAP, row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Clust4_ToClean) & Acinar_Seu_UMAP$UMAP_2 < 10 & Acinar_Seu_UMAP$UMAP_2 > -18 | row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Rest2) & Acinar_Seu_UMAP$UMAP_1 > -5 & Acinar_Seu_UMAP$UMAP_2 < 0)

Acinar_Clust5 = subset(Acinar_Seu_UMAP, row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Clust5_ToClean) & Acinar_Seu_UMAP$UMAP_2 > 0 | row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Rest2) & Acinar_Seu_UMAP$UMAP_1 > -15 & Acinar_Seu_UMAP$UMAP_2 > 10)


Acinar_Clust6 = subset(Acinar_Seu_UMAP, row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Clust6_ToClean) & Acinar_Seu_UMAP$UMAP_2 < -5  | row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Rest2) & Acinar_Seu_UMAP$UMAP_1 > -24 & Acinar_Seu_UMAP$UMAP_1 < -5 & Acinar_Seu_UMAP$UMAP_2 < -18)

Acinar_Clust8 = subset(Acinar_Seu_UMAP, row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Clust8_ToClean) & Acinar_Seu_UMAP$UMAP_2 < 5 & Acinar_Seu_UMAP$UMAP_1 < -10  | row.names(Acinar_Seu_UMAP) %in% colnames(Acinar_Rest2) &  Acinar_Seu_UMAP$UMAP_1 < -24 & Acinar_Seu_UMAP$UMAP_2 < -5)


Acinar_Clust7 = subset(Acinar_Seu, cells = c(row.names(Acinar_Clust1_Clean), row.names(Acinar_Clust2), row.names(Acinar_Clust3), row.names(Acinar_Clust4), row.names(Acinar_Clust5), row.names(Acinar_Clust6), row.names(Acinar_Clust8)), invert=T)

CheckInput = Acinar_Clust7
CheckUMAP(Acinar_Rest2)

set.res = c(1)
set.dim = c(10)
set.kparam = c(30)
ClusterFunc_All_RNA(Acinar_Rest2)



#Acinar_Clust6 = subset(Acinar_Seu, cells = c(row.names(Acinar_Clust1_Clean), row.names(Acinar_Clust2), row.names(Acinar_Clust3), row.names(Acinar_Clust4), row.names(Acinar_Clust5)), invert=T)


Acinar_Barcs_List = list("Acinar_1" = Acinar_Clust1_Clean, "Acinar_2" = Acinar_Clust2, "Acinar_3" = Acinar_Clust3, "Acinar_4" = Acinar_Clust4, "Acinar_5" = Acinar_Clust5, "Acinar_6" = Acinar_Clust6,  "Acinar_7" = Acinar_Clust7, "Acinar_8" = Acinar_Clust8)
Acinar_Barcs = GenerateMetaData(Acinar_Barcs_List)
###



### SUBCLUSTER DUCTAL #Use 2
Ductal_Seu = subset(Islet.combined, cells = row.names(Ductal_Clust))
DefaultAssay(Ductal_Seu) = "integrated"
Ductal_Seu <- RunPCA(Ductal_Seu, npcs = 10, verbose = FALSE)
Ductal_Seu <- RunUMAP(Ductal_Seu, reduction = "pca", dims = 1:10, spread=5)
Ductal_Seu_UMAP = as.data.frame(Ductal_Seu@reductions$umap@cell.embeddings)

set.res = c(1)
set.dim = c(10)
set.kparam = c(20)
ClusterFunc_All_RNA(Ductal_Seu)

DefaultAssay(Ductal_Seu) = "integrated"
Ductal_Seu <- FindNeighbors(Ductal_Seu, k.param=20, dims=1:10)
Ductal_Seu <- FindClusters(Ductal_Seu, resolution = 1)
Ductal_Clust1_ToClean = subset(Ductal_Seu, idents = c(4))
Ductal_Clust2_ToClean = subset(Ductal_Seu, idents = c(0,5))
Ductal_Clust3_ToClean = subset(Ductal_Seu, idents = c(3,7))
Ductal_Clust4_ToClean = subset(Ductal_Seu, idents = c(6))


Ductal_Clust1 = subset(Ductal_Seu_UMAP, row.names(Ductal_Seu_UMAP) %in% colnames(Ductal_Clust1_ToClean)  & Ductal_Seu_UMAP$UMAP_2 > 5 &  Ductal_Seu_UMAP$UMAP_1 < 7  | Ductal_Seu_UMAP$UMAP_1 > -9 & Ductal_Seu_UMAP$UMAP_1 < 6 & Ductal_Seu_UMAP$UMAP_2 > 6 )

Ductal_Clust2 = subset(Ductal_Seu_UMAP, row.names(Ductal_Seu_UMAP) %in% colnames(Ductal_Clust2_ToClean) & ! row.names(Ductal_Seu_UMAP) %in% row.names(Ductal_Clust1)  & Ductal_Seu_UMAP$UMAP_2 > -10 &  Ductal_Seu_UMAP$UMAP_1 > 3  | ! row.names(Ductal_Seu_UMAP) %in% row.names(Ductal_Clust1) & Ductal_Seu_UMAP$UMAP_1 > 10  & Ductal_Seu_UMAP$UMAP_2 > -9 )

Ductal_Clust3 = subset(Ductal_Seu_UMAP, row.names(Ductal_Seu_UMAP) %in% colnames(Ductal_Clust3_ToClean)  & Ductal_Seu_UMAP$UMAP_2 > -8 &  Ductal_Seu_UMAP$UMAP_1 < -10 )

Ductal_Clust4 = subset(Ductal_Seu_UMAP, row.names(Ductal_Seu_UMAP) %in% colnames(Ductal_Clust4_ToClean)  & Ductal_Seu_UMAP$UMAP_2 < -6 &  Ductal_Seu_UMAP$UMAP_1 > 3  | Ductal_Seu_UMAP$UMAP_2 < -9 &  Ductal_Seu_UMAP$UMAP_1 > 6 )

Ductal_Clust5 = subset(Ductal_Seu, cells = c(row.names(Ductal_Clust1), row.names(Ductal_Clust2), row.names(Ductal_Clust3), row.names(Ductal_Clust4)), invert=T)


CheckInput = Ductal_Clust5
CheckUMAP(Ductal_Seu)

Ductal_Barcs_List = list("Ductal_1" = Ductal_Clust1, "Ductal_2" = Ductal_Clust2, "Ductal_3" = Ductal_Clust3, "Ductal_4" = Ductal_Clust4, "Ductal_5" = Ductal_Clust5)
Ductal_Barcs = GenerateMetaData(Ductal_Barcs_List)
###

### SUBCLUSTER DELTA
Delta_Seu = subset(Islet.combined, cells = row.names(Delta_Clust))
DefaultAssay(Delta_Seu) = "integrated"
Delta_Seu <- RunPCA(Delta_Seu, npcs = 30, verbose = FALSE)
Delta_Seu <- RunUMAP(Delta_Seu, reduction = "pca", dims = 1:10, spread=5)
Delta_Seu_UMAP = as.data.frame(Delta_Seu@reductions$umap@cell.embeddings)

set.res = c(1)
set.dim = c(10)
set.kparam = c(50)
#ClusterFunc_All_RNA(Delta_Seu)

DefaultAssay(Delta_Seu) = "integrated"
Delta_Seu <- FindNeighbors(Delta_Seu, k.param=50, dims=1:10)
Delta_Seu <- FindClusters(Delta_Seu, resolution = 1)
Delta_Clust1_ToClean = subset(Delta_Seu, idents = c(2,3))

Delta_Clust1 = subset(Delta_Seu_UMAP, row.names(Delta_Seu_UMAP) %in% colnames(Delta_Clust1_ToClean) &  Delta_Seu_UMAP$UMAP_1 > 3 | Delta_Seu_UMAP$UMAP_1 > 8 )

CheckInput = Delta_Clust1
CheckUMAP(Delta_Seu)

Delta_Rest2 = subset(Delta_Seu, cells = row.names(Delta_Clust1), invert=T)

set.res = c(1)
set.dim = c(10)
set.kparam = c(50)
ClusterFunc_All_RNA(Delta_Rest2)

DefaultAssay(Delta_Rest2) = "integrated"
Delta_Rest2 <- FindNeighbors(Delta_Rest2, k.param=50, dims=1:10)
Delta_Rest2 <- FindClusters(Delta_Rest2, resolution = 1)
Delta_Clust2_ToClean = subset(Delta_Rest2, idents = c(2))
Delta_Clust3_ToClean = subset(Delta_Rest2, idents = c(3))

CheckInput = Delta_Clust2_ToClean
CheckUMAP(Delta_Rest2)

Delta_Clust2 = subset(Delta_Seu_UMAP, row.names(Delta_Seu_UMAP) %in% colnames(Delta_Clust2_ToClean) &  Delta_Seu_UMAP$UMAP_1 < -5 &  Delta_Seu_UMAP$UMAP_2 < 5 &  Delta_Seu_UMAP$UMAP_2 > -5 | row.names(Delta_Seu_UMAP) %in% colnames(Delta_Rest2) &  Delta_Seu_UMAP$UMAP_1 < -10 &  Delta_Seu_UMAP$UMAP_2 < 5 &  Delta_Seu_UMAP$UMAP_2 > -4)

Delta_Clust3 = subset(Delta_Seu_UMAP, row.names(Delta_Seu_UMAP) %in% colnames(Delta_Clust3_ToClean) & Delta_Seu_UMAP$UMAP_2 > 5 & Delta_Seu_UMAP$UMAP_1 < 2 | row.names(Delta_Seu_UMAP) %in% colnames(Delta_Rest2) &  Delta_Seu_UMAP$UMAP_2 > 5 & Delta_Seu_UMAP$UMAP_1 < -4 | row.names(Delta_Seu_UMAP) %in% colnames(Delta_Rest2) &  Delta_Seu_UMAP$UMAP_2 > 6.5 & Delta_Seu_UMAP$UMAP_1 < 0)



CheckInput = Delta_Clust3
CheckUMAP(Delta_Rest2)

Delta_Clust4 = subset(Delta_Seu, cells = c(row.names(Delta_Clust1), row.names(Delta_Clust2), row.names(Delta_Clust3)), invert=T)
CheckInput = Delta_Clust3
CheckUMAP(Delta_Seu)

Delta_Barcs_List = list("Delta_1" = Delta_Clust1, "Delta_2" = Delta_Clust2, "Delta_3" = Delta_Clust3, "Delta_4" = Delta_Clust4)
Delta_Barcs = GenerateMetaData(Delta_Barcs_List)
###


### SUBCLUSTER STELLATE #Use 24
Stellate_Seu = subset(Islet.combined, cells = row.names(Stellate_Clust))
DefaultAssay(Stellate_Seu) = "integrated"
Stellate_Seu <- RunPCA(Stellate_Seu, npcs = 40, verbose = FALSE)
Stellate_Seu <- RunUMAP(Stellate_Seu, reduction = "pca", dims = 1:20, spread=10)

Stellate_Barcs_List = list("Stellate" = Stellate_Seu)
Stellate_Barcs = GenerateMetaData(Stellate_Barcs_List)
###


### SUBCLUSTER ENDO
Endo_Seu = subset(Islet.combined, cells = row.names(Endo_Clust))
DefaultAssay(Endo_Seu) = "integrated"
Endo_Seu <- RunPCA(Endo_Seu, npcs = 30, verbose = FALSE)
Endo_Seu <- RunUMAP(Endo_Seu, reduction = "pca", dims = 1:10, spread=2)
Endo_Seu_UMAP = as.data.frame(Endo_Seu@reductions$umap@cell.embeddings)

set.res = c(1)
set.dim = c(10)
set.kparam = c(2)
ClusterFunc_All_RNA(Endo_Seu)


DefaultAssay(Endo_Seu) = "integrated"
Endo_Seu <- FindNeighbors(Endo_Seu, k.param=2, dims=1:10)
Endo_Seu <- FindClusters(Endo_Seu, resolution = 1)
ArterialCells = subset(Endo_Seu, idents = c(1,6,8,11))


ArterialCells2 = subset(Endo_Seu_UMAP, row.names(Endo_Seu_UMAP) %in% colnames(ArterialCells) | Endo_Seu_UMAP$UMAP_1 > 1.5 & Endo_Seu_UMAP$UMAP_2 < 2.5)
VenousCells = subset(Endo_Seu, cells = row.names(ArterialCells2), invert=T)

CheckInput = ArterialCells
CheckUMAP(Endo_Seu)

Endo_List = list("Arterial" = ArterialCells2, "Venous" = VenousCells)
Endo_Barcs = GenerateMetaData(Endo_List)

#No Further Subclusters - Only FLT1 + VWF cells
###

### SUBCLUSTER WBCs
WBCs_Seu = subset(Islet.combined, cells = row.names(WBC_Clust))
DefaultAssay(WBCs_Seu) = "integrated"
WBCs_Seu <- RunPCA(WBCs_Seu, npcs = 10, verbose = FALSE)
WBCs_Seu <- RunUMAP(WBCs_Seu, reduction = "pca", dims = 1:10, spread=10)
WBCs_Seu_UMAP = as.data.frame(WBCs_Seu@reductions$umap@cell.embeddings)

set.res = c(1)
set.dim = c(10)
set.kparam = c(10)
ClusterFunc_All_RNA(WBCs_Seu)

DefaultAssay(WBCs_Seu) = "integrated"
WBCs_Seu <- FindNeighbors(WBCs_Seu, k.param=10, dims=1:10)
WBCs_Seu <- FindClusters(WBCs_Seu, resolution = 1)
Mast_Cells = subset(WBCs_Seu, idents = c(2,4))
TCells = subset(WBCs_Seu, idents = c(0,1))
APCs = subset(WBCs_Seu, idents = c(3))


WBCs_Barcs_List = list("Mast_Cells" = Mast_Cells, "T_Cells" = TCells, "APCs" = APCs)
WBCs_Barcs = GenerateMetaData(WBCs_Barcs_List)
###


### SUBCLUSTER Gamma
Gamma_Seu = subset(Islet.combined, cells = row.names(Gamma_Clust))
DefaultAssay(Gamma_Seu) = "integrated"
Gamma_Seu <- RunPCA(Gamma_Seu, npcs = 30, verbose = FALSE)
Gamma_Seu <- RunUMAP(Gamma_Seu, reduction = "pca", dims = 1:20, spread=10)

Gamma_Barcs_List = list("Gamma" = Gamma_Seu)
Gamma_Barcs = GenerateMetaData(Gamma_Barcs_List)
###


### SUBCLUSTER Cycling 1
Cycling_Seu = subset(Islet.combined, cells = row.names(Cycling_Clust))
DefaultAssay(Cycling_Seu) = "integrated"
Cycling_Seu <- RunPCA(Cycling_Seu, npcs = 10, verbose = FALSE)
Cycling_Seu <- RunUMAP(Cycling_Seu, reduction = "pca", dims = 1:10, spread =  2)
Cycling_Seu_UMAP = as.data.frame(Cycling_Seu@reductions$umap@cell.embeddings)

set.res = c(1)
set.dim = c(10)
set.kparam = c(10)
ClusterFunc_All_RNA(Cycling_Seu)

DefaultAssay(Cycling_Seu) = "integrated"
Cycling_Seu <- FindNeighbors(Cycling_Seu, k.param=10, dims=1:10)
Cycling_Seu <- FindClusters(Cycling_Seu, resolution = 1)
Cycling_Clust_1_ToClean = subset(Cycling_Seu, idents = c(0))
Cycling_Clust_2_ToClean = subset(Cycling_Seu, idents = c(1))


Cycling_Clust1 = subset(Cycling_Seu_UMAP, row.names(Cycling_Seu_UMAP) %in% colnames(Cycling_Clust_1_ToClean) & Cycling_Seu_UMAP$UMAP_1 < -1.5 | Cycling_Seu_UMAP$UMAP_1 < 0 & Cycling_Seu_UMAP$UMAP_2 > -1 )


Cycling_Clust2 = subset(Cycling_Seu_UMAP, row.names(Cycling_Seu_UMAP) %in% colnames(Cycling_Clust_2_ToClean) & Cycling_Seu_UMAP$UMAP_2 > -1  | Cycling_Seu_UMAP$UMAP_2 > -0.5 & Cycling_Seu_UMAP$UMAP_1 > 0.5 & Cycling_Seu_UMAP$UMAP_1 < 7)
CheckInput = Cycling_Clust2
CheckUMAP(Cycling_Seu)

Cycling_Clust3 = subset(Cycling_Seu, cells = c(row.names(Cycling_Clust1), row.names(Cycling_Clust2)), invert=T)

Cycling_Barcs_List = list("Cycling_1" = Cycling_Clust1, "Cycling_2" = Cycling_Clust2, "Cycling_3" = Cycling_Clust3)
Cycling_Barcs = GenerateMetaData(Cycling_Barcs_List)
###


### SUBCLUSTER Macro
Macro_Seu = subset(Islet.combined, cells = row.names(Macro_Clust))
DefaultAssay(Macro_Seu) = "integrated"
Macro_Seu <- RunPCA(Macro_Seu, npcs = 30, verbose = FALSE)
Macro_Seu <- RunUMAP(Macro_Seu, reduction = "pca", dims = 1:10, spread=5)

Macro_Barcs_List = list("Macro" = Macro_Seu)
Macro_Barcs = GenerateMetaData(Macro_Barcs_List)
###

### SUBCLUSTER Epsilon
Epsilon_Seu = subset(Islet.combined, cells = row.names(Epsilon_Clust))
DefaultAssay(Epsilon_Seu) = "integrated"
Epsilon_Seu <- RunPCA(Epsilon_Seu, npcs = 10, verbose = FALSE,n.neighbors=5)
Epsilon_Seu <- RunUMAP(Epsilon_Seu, reduction = "pca", dims = 1:10,verbose = F,n.neighbors=5)

Epsilon_Barcs_List = list("Epsilon" = Epsilon_Seu)
Epsilon_Barcs = GenerateMetaData(Epsilon_Barcs_List)

Generate plots

List_Meta = list("Alpha" = Alpha_Barcs, "Beta" = Beta_Barcs,"Delta" = Delta_Barcs,  "Epsilon" = Epsilon_Barcs, "Gamma" = Gamma_Barcs, "Acinar" = Acinar_Barcs, "Ductal" = Ductal_Barcs, "Stellate" = Stellate_Barcs, "Cycling" = Cycling_Barcs, "Macro" = Macro_Barcs,   "Endo" = Endo_Barcs, "WBCs" = WBCs_Barcs)

List_Seu = list( "Alpha" = Alpha_Seu, "Beta" = Beta_Seu, "Delta" =  Delta_Seu, "Epsilon" = Epsilon_Seu, "Gamma" = Gamma_Seu,"Acinar" = Acinar_Seu, "Ductal" =  Ductal_Seu, "Stellate" = Stellate_Seu, "Cycling" = Cycling_Seu,"Macrophages"= Macro_Seu, "Endothelium" = Endo_Seu, "WBCs" = WBCs_Seu)


sz = 8
DPList = list()

Islet.combined = AddMetaData(Islet.combined, Islets_Barcs, "Islets_Barcs")

Idents(Islet.combined) = "Islets_Barcs"
DPList[["AllCells"]] = DimPlot(Islet.combined, label = T, repel=T, label.size = sz*0.3) + NoLegend()+theme(axis.title = element_text(size = sz), axis.text = element_text(size = sz))

Barcodes_Granular = as.data.frame(matrix(ncol=1, nrow=0))
colnames(Barcodes_Granular) = "Pop"
sz = 8
for(x in 1:length(List_Seu)){
Pull = List_Seu[[x]]
PullMeta = List_Meta[[x]]
Pull = AddMetaData(Pull, List_Meta[[x]], "Assigns")
Idents(Pull) = "Assigns"
DPList[[names(List_Seu[x])]] = DimPlot(Pull, label = T, label.size = sz*0.3) + NoLegend() + theme(axis.title = element_text(size = sz*0.7), axis.text = element_text(size = sz*0.7), plot.title = element_text(size = sz, face="bold")) + ggtitle(names(List_Seu[x]))
Barcodes_Granular = rbind(Barcodes_Granular, PullMeta)
}
write.csv(Barcodes_Granular, "Islets_Barcodes_Granular_5AUG22.csv")

Genes = unique(c("GCG", "IRX1", "IRX2", "FAP", "DPP4", "LOXL4", "PLCE1", "NEUROD1", "PPP1R1A", "INS", "NKX6-1", "DLK1", "ADCYAP1", "ABCC8", "RBP4", "MAFA", "FFAR4", "PDX1","NKX6-1", "SIX3", "SST", "NPY", "POU3F1", "LEPR", "GHRL","ETV1",  "PPY", "PRSS1", "CPA1", "CPB1", "CTRC", "PLA2G1B","SDC4", "LCN2",  "KRT19", "MUC1", "AMBP", "COL1A1", "DKK3", "KCNE4", "COL15A1", "RGS5", "SPARC", "VCAN","PDGFRB","SERPINH1", "TBX3", "COL1A1", "SOX10", "MKI67", "UBE2C", "CDK1", "CENPA", "FAM111B",  "SDS", "MERTK", 'ITGAX', "VWF", "CLDN5", "FLT1", "ACTA2", "RGS5", "LUM", "DCN", "TPSAB1", "KIT", "MS4A2", "CD3E", "TRAC", "IL7R", "CCL17", "IDO1"))
write.csv(as.data.frame(table(Barcodes_Granular$Pop)), "Islet_Clusters_5AUG22.csv")


Barcodes_Granular$Pop = factor(Barcodes_Granular$Pop, levels = rev(c("Alpha_1", "Alpha_2", "Alpha_3", "Alpha_4","Alpha_5", "Alpha_6", "Beta_1", "Beta_2", "Beta_3", "Beta_4", "Beta_5", "Beta_6", "Beta_7", "Beta_8", "Beta_9", "Delta_1", "Delta_2", "Delta_3", "Delta_4", "Epsilon", "Gamma", "Acinar_1", "Acinar_2", "Acinar_3", "Acinar_4", "Acinar_5", "Acinar_6", "Acinar_7", "Acinar_8", "Ductal_1", "Ductal_2", "Ductal_3", "Ductal_4", "Ductal_5", "Stellate", "Cycling_1", "Cycling_2", "Cycling_3", "Macro", "Arterial","Venous", "Mast_Cells", "T_Cells", "APCs")))


Islet.combined = AddMetaData(Islet.combined, Barcodes_Granular, "Gran")
Idents(Islet.combined) = "Gran"
DefaultAssay(Islet.combined) = "RNA"
myPalette <- colorRampPalette(c("lightgrey", "navy"))

Dotz = DotPlot(object = Islet.combined, features = Genes, assay = "RNA",  cols = c("lightgrey", "navy"), dot.scale = 1.5)  + theme(legend.position = "bottom", axis.text.x = element_text(angle = 90, hjust = 1, vjust=0.5, size=sz*0.7, face = "italic"), axis.text.y = element_text(size=sz*0.7), legend.text = element_text(size=sz), legend.title = element_text(size=sz)) + xlab("") + ylab("")+ scale_colour_gradientn(colours = myPalette(100), limits = c(0,2.5), na.value = "white")
#
pdf(paste("ISLETS_DotPlot_RNA_5AUG22.pdf", sep=""), width=20, height=10)
print(Dotz)
dev.off()

Patch = DPList[[1]] + DPList[[2]] + DPList[3] + DPList[4] + DPList[5] + DPList[6] + DPList[7] + DPList[8] + DPList[9] + DPList[10] + DPList[11] + DPList[12] + DPList[13] + Dotz +
   plot_layout(design = "AAAABBCC
            AAAABBCC
            AAAADDEE
            AAAADDEE
            FFGGHHII
            FFGGHHII
            JJKKLLMM
            JJKKLLMM
            OOOOOOOO
            OOOOOOOO
            OOOOOOOO
            OOOOOOOO
            OOOOOOOO")
pdf("Islets_Figure_RNA_5AUG22.pdf", width = 8, height = 12)
print(Patch)
dev.off()


   plot_layout(design = "AABC
              AADE
              FGHI
              JKLM
             OOOO
             OOOO")

Islet.combinedData = as.data.frame(t(Islet.combined@assays$RNA@counts))
Islet.combinedData_subs = Islet.combinedData %>% dplyr::select("GCG", "IRX1", "INS", "DLK1", "SST", "POU3F1")

Islet.combinedDecontx = as.data.frame(t(Islet.combined@assays$decont@counts))
Islet.combinedDecontx_subs = Islet.combinedDecontx %>% dplyr::select("GCG", "IRX1", "INS", "DLK1", "SST", "POU3F1")

for(b in unique(Islets_Barcs$Pop)){
if(b == "Stellate"){PlotGenes = c("COL1A1", GeneLists[["Baron_Stellate"]])}else if(b == "Endo"){PlotGenes = c("CSPG4", "GJA4", "APLN", "CD82", "CHST1", "PTH1R", "BMX", "GKN3", "IGFBP2", "DCN", "LIM", "PDGFRA", "IL33", GeneLists[["EndoList"]])}else if(b == "Beta"){PlotGenes = c("INS", GeneLists[["Seger_Beta"]], GeneLists[["Beta_Subpops"]])}else if(b == "Acinar"){PlotGenes = c("PRSS1", GeneLists[["Seger_Acinar"]])}else if(b == "Ductal"){PlotGenes = c(GeneLists[["Peng_Ductal"]])}else if(b == "Alpha"){PlotGenes = c("GCG", GeneLists[["Seger_Alpha1"]], GeneLists[["Seger_Alpha2"]])}else if(b == "Macro"){PlotGenes = c(GeneLists[["MacroCells"]])}else if(b == "WBCs"){PlotGenes = c(GeneLists[["MainList"]], "TPSAB1", GeneLists[["MastCells"]])}else if(b == "Delta"){PlotGenes = c("SST", GeneLists[["Rorsman_Pardo_Delta"]])}else if(b == "Gamma"){PlotGenes = c("PPY" , GeneLists[["GammaCells"]])}

pdf(paste("Islets_DotPlots_", b, "MarkerGenes_3AUG22.pdf", sep=""), width = 12, height = length(PlotGenes)/4)
print(DotPlot(object = Islet.combined, features = unique(PlotGenes), assay = "RNA", cluster.idents = T)+coord_flip()+theme(axis.text.x = element_text(angle=90, hjust=1, vjust = 0.5)))
dev.off()
paste(c(b, ":", PlotGenes), collapse = ", ")
}




### Barcodes for GlowwormInput
Islets_Barcodes = Barcodes_Granular
Islets_Barcodes$Broad = gsub("T_Cells", "WBCs", gsub("APCs", "WBCs", gsub("Mast_Cells", "WBCs", gsub("Arterial","Endothelial", gsub("Venous", "Endothelial", Islets_Barcodes$Pop)))))
Islets_Barcodes$Broad = gsub("_.*", "",Islets_Barcodes$Broad)

save(list=c("Islets_Barcodes", "Islet.combined"), file = "~/Islets_Glowworm_15AUG22.RData")

#######



sz = 8
DPList = list()

Idents(Islet.combined) = "Islets_Barcs"
DPList[["AllCells"]] = DimPlot(Islet.combined) + NoLegend()+theme(axis.title = element_text(size = sz), axis.text = element_blank(), axis.ticks = element_blank())

sz = 8
for(x in 1:length(List_Seu)){
Pull = List_Seu[[x]]
PullMeta = List_Meta[[x]]
Pull = AddMetaData(Pull, List_Meta[[x]], "Assigns")
Idents(Pull) = "Assigns"
PullName = names(List_Seu[x])
DPList[[names(List_Seu[x])]] = DimPlot(Pull) + NoLegend() + theme(axis.title = element_blank(), axis.text = element_blank(), plot.title = element_text(size = sz, face="bold"), axis.ticks = element_blank()) + ggtitle(paste(PullName, " (",length(unique(PullMeta$Pop)), ")",  sep=""))
}

Patch = wrap_plots(DPList) +
   plot_layout(design = "AAABCDE
                         AAAFGHI
                         AAAJKLM")
pdf("Islets_Figure_RNA_SEP22.pdf", width = 20, height = 8)
print(Patch)
desv.off()

pdf("Fetal_Legends.pdf")
print(DimPlot(Fetal.combined, label=T))
dev.off()

List_Seu = list( "Alpha" = Alpha_Seu, "Beta" = Beta_Seu, "Delta" =  Delta_Seu, "Epsilon" = Epsilon_Seu, "Gamma" = Gamma_Seu,"Acinar" = Acinar_Seu, "Ductal" =  Ductal_Seu, "Stellate" = Stellate_Seu, "Cycling" = Cycling_Seu,"Macrophages"= Macro_Seu, "Endothelium" = Endo_Seu, "WBCs" = WBCs_Seu)


#Clean up Beta
DefaultAssay(Beta_Seu) = "RNA"
FP = FeaturePlot(Beta_Seu, unique(c(GeneLists[["Seger_Beta"]], GeneLists[["Beta_Subpops"]], GeneLists[["Lori_Extras"]])), ncol = 5)

pdf("Beta_Clust_FeaturePlots_MAR23.pdf", width = 14, height = 2*ceiling(length( unique(c(GeneLists[["Seger_Beta"]], GeneLists[["Beta_Subpops"]], GeneLists[["Lori_Extras"]])))/5))
print(FP)
dev.off()

Beta_Seu = AddMetaData(Beta_Seu, Beta_Barcs, "BETA")
Idents(Beta_Seu)= "BETA"
VlnPlot(Beta_Seu, c("INS", "FFAR4", "RBP4", "ID1", "ID2", "ID3"), pt.size = 0)

DimPlot(Beta_Seu, label = T)

Beta_Barcs_V2 = Beta_Barcs
Beta_Barcs_V2$Pop = gsub("Beta_7", "Beta_6", gsub("Beta_8", "Beta_5", gsub("Beta_6", "Beta_4", gsub("Beta_9", "Beta_4",  Beta_Barcs_V2$Pop))))



#Clean up Alpha
DefaultAssay(Alpha_Seu) = "RNA"
FP = FeaturePlot(Alpha_Seu, unique(c(GeneLists[["Seger_Alpha1"]], GeneLists[["Seger_Alpha2"]], GeneLists[["Lori_Extras"]])), ncol = 5)

pdf("Alpha_Clust_FeaturePlots_MAR23.pdf", width = 14, height = 2*ceiling(length( unique(c(GeneLists[["Seger_Alpha1"]], GeneLists[["Seger_Alpha2"]], GeneLists[["Lori_Extras"]])))/5))
print(FP)
dev.off()

Alpha_Seu = AddMetaData(Alpha_Seu, Alpha_Barcs, "Alpha")
Idents(Alpha_Seu)= "Alpha"

DimPlot(Alpha_Seu, label = T)

Alpha_Barcs_V2 = Alpha_Barcs
Alpha_Barcs_V2$Pop = gsub("Alpha_4", "Alpha_3", gsub("Alpha_5", "Alpha_3", gsub("Alpha_6", "Alpha_3",   Alpha_Barcs_V2$Pop)))


#Clean up Delta
DefaultAssay(Delta_Seu) = "RNA"
FP = FeaturePlot(Delta_Seu, unique(c(GeneLists[["Rorsman_Pardo_Delta"]],  GeneLists[["Lori_Extras"]])), ncol = 5)

pdf("Delta_Clust_FeaturePlots_MAR23.pdf", width = 14, height = 2*ceiling(length( unique(c(GeneLists[["Rorsman_Pardo_Delta"]], GeneLists[["Lori_Extras"]])))/5))
print(FP)
dev.off()

Delta_Seu = AddMetaData(Delta_Seu, Delta_Barcs, "Delta")
Idents(Delta_Seu)= "Delta"

DimPlot(Delta_Seu, label = T)

Delta_Barcs_V2 = Delta_Barcs
Delta_Barcs_V2$Pop = gsub("Delta_4", "Delta_1",  Delta_Barcs_V2$Pop)


#Clean up Acinar
DefaultAssay(Acinar_Seu) = "RNA"
FP = FeaturePlot(Acinar_Seu, unique(c(GeneLists[["Seger_Acinar"]], GeneLists[["Lori_Extras"]])), ncol = 5)

pdf("Acinar_Clust_FeaturePlots_MAR23.pdf", width = 14, height = 2*ceiling(length( unique(c(GeneLists[["Seger_Acinar"]],  GeneLists[["Lori_Extras"]])))/5))
print(FP)
dev.off()

Acinar_Seu = AddMetaData(Acinar_Seu, Acinar_Barcs, "Acinar")
Idents(Acinar_Seu)= "Acinar"

DimPlot(Acinar_Seu, label = T)


Acinar_Reclust = subset(Acinar_Seu, idents = c("Acinar_6", "Acinar_4"))

DefaultAssay(Acinar_Reclust) = "integrated"
Acinar_Reclust <- FindNeighbors(Acinar_Reclust, k.param=30, dims=1:10)
Acinar_Reclust <- FindClusters(Acinar_Reclust, resolution = 1)
DimPlot(Acinar_Reclust, label = T)
FeaturePlot(Acinar_Seu, "KRT19")
Acinar_Cluster3 = subset(Acinar_Reclust, idents = c(0,2,5))


Acinar_Barcs_V2 = Acinar_Barcs
Acinar_Barcs_V2$Pop = gsub("Acinar_3", "Acinar_2", gsub("Acinar_5", "Acinar_2",gsub("Acinar_4", "Acinar_2", gsub("Acinar_6", "Acinar_2",  Acinar_Barcs_V2$Pop))))
Acinar_Barcs_V2$Pop = ifelse(row.names(Acinar_Barcs_V2) %in% colnames(Acinar_Cluster3), "Acinar_3", Acinar_Barcs_V2$Pop)

Acinar_Seu = AddMetaData(Acinar_Seu, Acinar_Barcs_V2, "Acinar")
Idents(Acinar_Seu)= "Acinar"

DimPlot(Acinar_Seu, label = T)

#Clean up Ductal
DefaultAssay(Ductal_Seu) = "RNA"
FP = FeaturePlot(Ductal_Seu, unique(c(GeneLists[["Peng_Ductal"]],  GeneLists[["Lori_Extras"]])), ncol = 5)

pdf("Ductal_Clust_FeaturePlots_MAR23.pdf", width = 14, height = 2*ceiling(length( unique(c(GeneLists[["Seger_Ductal1"]], GeneLists[["Lori_Extras"]])))/5))
print(FP)
dev.off()


FeaturePlot(Ductal_Seu, c("AMBP","FXYD2", "MUC1", "FXYD3",  "KRT19"))

Ductal_Seu = AddMetaData(Ductal_Seu, Ductal_Barcs, "Ductal")
Idents(Ductal_Seu)= "Ductal"

DimPlot(Ductal_Seu, label = T)

Ductal_Barcs_V2 = Ductal_Barcs
Ductal_Barcs_V2$Pop = gsub("Ductal_4", "Ductal_3", gsub("Ductal_5", "Ductal_3",  Ductal_Barcs_V2$Pop))
###


Cycling_Barcs_V2 = Cycling_Barcs
Cycling_Barcs_V2$Pop = gsub("Cycling_1", "Beta_7 [Cycling]", gsub("Cycling_2", "Beta_8 [Cycling]", gsub("Cycling_3", "Alpha_4 [Cycling]", Cycling_Barcs_V2$Pop)))


Islet_Assigns_MAR23 = rbind(Ductal_Barcs_V2, Delta_Barcs_V2, Acinar_Barcs_V2, Alpha_Barcs_V2, Beta_Barcs_V2, Epsilon_Barcs, Gamma_Barcs, Stellate_Barcs, Cycling_Barcs_V2, Macro_Barcs,  Endo_Barcs, WBCs_Barcs)

Islet_Assigns_MAR23$Broad = Islet_Assigns_MAR23$Pop
Islet_Assigns_MAR23$Broad = gsub("T_Cells", "WBCs", gsub("APCs", "WBCs", gsub("Mast_Cells", "WBCs", gsub("Arterial","Endothelial", gsub("Venous", "Endothelial", Islet_Assigns_MAR23$Pop)))))
Islet_Assigns_MAR23$Broad = gsub("_.*", "",Islet_Assigns_MAR23$Broad)
write.csv(Islet_Assigns_MAR23, "Islet_Assigns_MAR23.csv")

save(list=c("Islet_Assigns_MAR23", "Islet.combined"), file = "~/Islets_Glowworm_28MAR23.RData")


Hannahglover/Glowworm documentation built on Jan. 16, 2024, 11:47 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close